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ABSTRACT 

We inti-oduce the star cluster evolution code EVOLVE Me A CLUSTER OF StarS 
(EMACSS), a simple yet physically motivated computational model that describes the evo- 
lution of some fundamental properties of star clusters in static tidal fields. We base our pre- 
scription upon the flow of energy within the cluster, which is a constant fraction of the total 
energy per half-mass relaxation time. According to Henon's predictions, this flow is indepen- 
dent of the precise mechanisms for energy production within the core, and therefore does not 
require a complete description of the many-body interactions therein. For a cluster of equal- 
mass stars, we thence use dynamical theory and analytic descriptions of escape mechanisms 
to construct a series of coupled differential equations expressing the time-evolution of cluster 
mass and radius. These equations are numerically solved using a 4th order Runge-Kutta inte- 
gration kernel, and the results bench-marked against a database of direct A^-body simulations. 
We use simulations containing a modest initial number of stars (1024 ^ A^ ^ 65536), and 
point-mass tidal fields of various strengths. Our prescription is publicly available, and repro- 
duces the A^-body results to within ^ 10% accuracy for the entire post-collapse evolution of 
star clusters. 

Key words: stellar dynamics: methods - galaxies: star clusters - globular clusters: general - 
methods: iV-body simulations -methods: Numerical 



1 INTRODUCTION 

The evol ution of a star cluster is driven by a comb i nation of re - 
laxation JAmbartsumianll 19381; Ichandrasekhaij|l942l; lKing|ll958l) . 



, binary in teractions dHeg gidl 1 9751) , stellar evolution and stellar en- 
■ counters JHut et al J 19921) . Furthermore, for a cluster located within 
the tidal field of a galaxy, effects resul ting from clu s ter's interaction 
with the tidal field can be important JHenonll 19601 ; iLee & Ostriken 
Il987h . with the conseq uences of a limiting 'Jacobi' radius, and a 
shortened total lifetime dBaumgardt & Makindl2003n . The result is 
a complex system, in which the simultaneous modelling of several 
properties (e.g. mass, half-mass radius, density profile) is equivo- 
cal. 

As a result of this complexity, dynamical simulations are an 
appealing manner through which to study the evolution of star clus- 
ters. Significant success has been achieved through direct TV-body 
integrations (see Aarseth 1973; Makino 1996; Spurzem 1999), al- 
beit at the expense of high co mputational cost. Likewise, faster 
schemes (Monte Carlo schemes; HenonI 19751 ; lGierszlll998|, solv- 
ing th e Fokker-Planck equation; ICohnlll979l Gas models; iLarsonl 
Il970h have each obtained significant success modelling high A'^- 
systems, although with a reduction in versatility owing to the as- 
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sumptions required for each (IHeggie & Hutil2003h . Of particular 
relevance to this study, the Fokker-Planck equation represents one 
successful and explicit formula for modelling the evolution of the 
distribution function of a many-body system, in a manner similar 
to our intended prescription. However, this formula generally does 
not unambiguously exp ress the dynamical properties for wh ich we 
intend our prescription dTakahashi & Portegies ZwardEOOOl) . 

Qualitatively, the long-term evolution of a star cluster has been 
characterised into three phases; core collapse, core bounce and ex- 
pansion, and tidally limited contraction. The first stage (core col- 
lapse) has historically been the most studied, and results from the 
diffusion of kinetic energy from core stars outward through re- 
laxation. These stars, upon loosing kinetic energy, move inward 
and thus experience a deeper potential. As a result, these inward- 
moving stars accelerate, eventually leading to a net increase in core 
kinetic energy (a consequence encapsulated by the concept of neg- 
ative heat capacity in gravothermal systems, iLvnden-Bell & Wood 
|l968) . However, the increasing kinetic energy of core stars will ul- 
timately enhance the rate at which relaxation diffuses energy from 
the core, which in turn accelerates the process of collapse. 

Core collaps e is eventually halted once a source of energ y be- 
comes viable (see lStatleret al.ll 19871; lOiersz & Heggielll994l) . and 
energy removed from the core can be replaced without the need for 
further contraction. At this point, the core ' bounces' outward owing 
to excess energy released during collapse dlnagaki & Lvnden-Belll 
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1 1983b . before (for an isolated or compact cluster) entering an ex- 
pansion ph ase wherein outflowing energy from the core inflates 
the cluster ( lHenoij|l96^ : lLightman & Shapirdll978l) . During this 
process, the rate of energy production in the core comes into 
balance with the flow of energy that is required for the global 
evolution. A final contraction stage occurs once the system has 
expanded sufficiently that the evaporation of stars over the Ja- 
cobi radius has become the dominant evolutionary consequence 
JGieles & Baumgardt 2008) . and the evolution of the system is de- 
fine d by the rate at whic h stars are lost to the tidal field. 

iGieles et al.l ( 1201 ih have proposed a formalism through which 
the expansion and tidally limited phases can be linked, and thus 
a complete description of the life cycle obtained. The two phases 
therefore constitute extreme cases, defining asymptotic behaviour 
in the evo lution of clust ers. Making the assumption of self-similar 
evolution. iHenonl (11965) showed that, for expanding clusters (with- 
out escaping stars) the half-mass relaxation time scales linearly 
with time (t,h oc t), and therefo re the h a lf-ma ss radius scales with 
t^'"^. By similar methodology, iHenonl ( 119611) predicted that if a 
cluster is limited by it's Jacobi radius, the mass decreases linearly in 
time, and the half-mass radius will scale linearly with the Jacobi ra- 
dius. Hence, the half-mass radius will scale with (icv — i)^' '^ where 
fev is the total lifetime, or evaporation time of a cluster These two 
power-laws constitute extremes of the evolution, with much of a 
cluster's life cycle forming a transition between these two phases. 

For the two interpretations quoted above, Henon considered 
evolut ion to occur in a self-similar fashion (described in Henon 
( 1196 ll) as being homologous); the shape of the density profile re- 
mains constant, with an isothermal cusp and finite truncation at the 
Jacobi radius. Accordingly, if a cluster evolves self-similarly, the 
dynamical effects are limited to a variation in overall scale, and a 
corresponding gradual reduction of total mass. Although this sim- 
plifying assumption is not valid in extreme cases (prior to the es- 
tablishment of balanced evolution, or when close to final dissolu- 
tion and containing a relatively low number of stars), we initially 
employ this interpretation before exploring the effects of none self- 
similar evolution. 

The goal of this paper is to present the initial stage in the devel- 
opment of a dynamical prescription encompassing both the expan- 
sion and tidal contraction phases. To this end, our initial prescrip- 
tion is strongly simplified; we consider only the dynamical effects 
stemming from the interaction of equal-mass stars without inter- 
nal evolution, and thereupon eliminate the complications of mass 
segregation. Although these effects are important for realistic glob- 
ular clusters (Spitzer 1969, 1975), such effects serve here mainly 
to complicate the dynamics we seek to study, and anyway do not 
necessarily i ntroduce effects not otherwise r epresented in an equal 
mass model dLvnden-Bell & Eggletonll98d) . We allow our models 
no primordial binary content, and model our stars as point particles 
so as to eliminate perturbations caused by direct stellar collisions. 
Our approximation for the tidal field into which we immerse the 
cluster is also simplified, as we model the tidal field to be that of a 
point mass galaxy located well outside the cluster The ensuing Ja- 
cobi radius is regarded as a spherical surface, although we do take 
into account the effect of preferential trajectories for escape (see 
section |23t . The model is constructed in the form of an efficient 
C++ codqj, which we calibrate against A'^-body simulations. 

In section|2]we discuss the mechanics and physical details of 
escape considered by our prescription, which is itself discussed in 
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section[3] Following this, in section|4]we discuss our TV-body sim- 
ulations, against which our prescription is calibrated and verified in 
section|5] Finally, we consider the success and shortcomings of our 
prescription in section |6] and outline the future physics that is to be 
incorporated into our model. 



2 EVOLUTION OF STAR CLUSTERS 

The dynamical evolution of a star cluster is a process primarily 
drive n by the radial diffusion of energy from the ' hot' (energetic) 
core dvon Hoemen 1957c IHenonl 1 96 ll : iLarsorJ 1 9701) . Thus, the sys- 
tem's life cycle is spent striving to establish equipartition. Equipar- 
tition remains elusive however, on account of the inherent nega- 
tive heat capacity of gravothermal systems, with the result that the 
flow of energy continually brings the cluster further out of equilib- 
rium. This means that the energy 'produced' in the corqjis contin- 
ually transferred by relaxation outward into the halo of the cluster, 
where the injection of energy drives the dynamical effects we seek 
to study. 

We shall examine star clusters both within a tidal field, and 
isolated from external tidal influences. Although such isolated clus- 
ters do not naturally exist, models of isolated clusters describe the 
expansion phase experi enced in the e arly evolution of (realistic) 
tidally limited cluster^ JAarsetlJll97ll) . 

We begin by taking the familiar expression for the virial radius 
of a star cluster, r = —GhP/2W, where W is potential energy, G 
is the gravitational constant and M the total mass. Assuming that 
the system remains in virial equilibrium throughout its evolution 
such that the total energy E — W/2, we can express this total 
energy as 

GN'^fh^ 



E 



4r 



(1) 



where Ai = Nfh, with A'^ the number of stars and fh the mean 
mass of stars. Alternatively, we can express energy in terms of the 
half-mass radius ri,, such that 



E = -K- 



GN^m^ 



Th 



(2) 



in which k is a form factor dependent upon the density profile 
whereby rh — inr. However, here we shall make no distinction, 
i.e. we assume k — 1/4. 

From these two expressions, we can begin to parametrise the 
implications of the outward diffusion of energy by making two ini- 
tial assumptions. Firstly, we assume that the flux of energy passing 
through any shell is equal to the change in energy inside that shell. 
Hence, if a cluster has achieved balanced evolution, the flux of en- 
ergy passing through a spherical shell located at any given radius 
is equivalent to that released form the core. Secondly, we assume 
that the fraction of the total energy passing through this shell is 
constant per relaxation time, as the system is restored to energetic 
equilibrium over this timescale. Thus we find 



E 



1 



(3) 



^ Principally by dynamical mechanisms suc h as the forma tion and harden- 
ing of binaiies during three body encounters ('H eggielll975b. a lthough other 
options such as st ellar evolution are v iable (Giel es et al.l2010l) . 
•^ To quote from Aarseth & Heggij 1 11998) . "much may be learned by the 
study of more tractable, idealised models, provided that the goal of under- 
standing the behaviour of real clusters is always kept in mind." 
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where ^ is a dimensionless constant and trh is the half-mass re- 
laxation time. In ac cordance with the derivation presented by 
ISpitzer & HartI ( 1197 ih . we define trh as 



t,-h = 0.138- 



/mG ln(7iV)' 



(4) 



in which ln(77V) is the C oulomb logarithm wit h 7 ~ 0.11 
for equal mass clusters fsee lGiersz & Heggidll994t) . Using equa- 
tion (TJ and noting that m is constant for an equal-mass system, it 
is evident that any change in energy will depend only upon A^ and 
r. Hence 



E_ _ N_ r_ 
LEI "" iV"^r' 



(5) 



while equation l|2j would result in an equivalent expression in 
which k/ hi would be included. However, as we have fixed k = 1/4, 
we approximate the evolution of both virial and half- mass radius by 
equation l|5j- 

Equation ^ demonstrates that a change of energy can have 
two principle dynamical effects; either stars are lost, carrying en- 
ergy away out of system, and/or the radius will vary. It is appro- 
priate at this stage is to parametrise the change in N and r per 
trh in terms of dimensionless escape and expansion rates ^ and 
/i ( lGoodmanlll984 ISaumgardt et al.ll2002h . From the definitions 
therein we therefore take. 
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If only one phase of evolution is modeled, both ^ and /j, can be 
considered to be constant (ibid.). By contrast, since we attempt to 
model both expansion and (in a tidal field) contraction concurrently, 
we must allow ^ and /i to vary throughout the life cycle. By com- 
bining equations ([3), ((SJl and Q, we find (^ is related to ^ and /i 
by 



C = M + 2e, 



(8) 



at all times throughout the life cycle. We hence show that the time- 
evolution of A'' and r will be defined by a set of differential equa- 
tions expressing E, N and f in terms of (", A^ and r. We now exam- 
ine these equations for a cluster in isolation (thereupon only experi- 
encing it's expansion phase, section |2T| ), or a tidal field (and hence 
experiencing both expansion and contraction, sections lZ2l and [23l >. 
and combine these in section[3l 



2.1 Isolated clusters 

Both expansion and mass-loss occur in isolated clusters 
(Baumgardt et al.l2002h meaning that equations ^ and l[7) must be 
solved simultaneously. Dividing equation ([TJ by equation JSJ and 
eliminating /i with equation {8}, we obtain the differential equa- 
tion. 



(9) 



dN N \ ^ 

It has been shown JBaumgardt et alj I2OO2I and references con- 
tained) that an isolated cluster will expand nearly (although not 
entirely) homologously throughout it's lifetime, loosing a roughly 
constant fraction of it's stars per relaxation time. As such, it is 



permissible to attempt a solution with constant 5. Making this as- 
sumption and writing the constant ^ of isolated clusters as ^1 , equa- 
tion {9)1 is separable and (following integration) yields 



r 



N_ 



_c_ 

5l 



(10) 



where TVo and ro are scale constants, usually considered to corre- 
spond to the cluster at the time of core collapse. As we have as- 
sumed r — Th, we can eliminate the rh dependence in relaxation 
time such that 



trh = trh.O 

where. 



- / ln(7Afo) \ 



2Ci 



3C " 7Ci 



(11) 



(12) 



If we now make the assumption that the Coulomb logarithm is con- 
stant, we can explicitly compute the evolution of A^ and r by ap- 
plying equation jilt to l[6j and solving. Thus, 






1 + 
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while, by a similar method, 

6 



-= 1 + 



i^trh.O 



(t - to) 



(13) 



(14) 



Equations ( |13t and ( 114b demonstrate a power-law scaling, and are 
dependent upon the time at which evolution begins, to. Choosing 
this time as. 



to = — trh.O, 

we obtain. 






ro 



^ s (2+i.)/3 



(15) 



(16) 



(17) 



which are identical to the forms suggested by lOoodmanl (119841) . 
Our chosen definition of to, equation dlSt is here used as a numer- 
ical tool to simplify equations l lI3t and ( I14t . although implicitly 
assumes that the solution of A'^ and r will pass through the origin at 
t = 0. This point (to) corresponds to the moment at which a cluster 
enters balanced evolution, and roughly to the core collapse time of 
the system. 

Equations ( I16t and ( I17l l imply that the curves A'^(t) and r(f) 
are parallel in a log-log figure for a cluster of any A^o (see fig- 
ure[TJa)). Inclusion of the A'^ dependence of the Coulomb logarithm 
(numerically, since the logarithmic terms retained in equation dlSt 
prohibit integration of equation (|7j) modifies this evolution to al- 
low a limited convergence of the evolutionary tracks, which is es- 
pecially appar ent at late times (fi gure [jib)). 

However. iBaumgardt et al.M 2002) show (via A'-body simula- 
tions) that there is in fact a crossing of the evolutionary tracks (i.e. 
^ = S,{N), with larger ^ for larger A'^). This effect is most likely 
a result of a positive correlation b etween core density and A'^ for a 
star cluster dHeggie & Hulll2003r) . Stars escaping from an isolated 
cluster will by necessity escape to infinity, as there is no possibil- 
ity of escape via translation over a tidal boundary. It follows that 
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10° 10^ 10" 10" lO" 10" 10" 10" 10" 10" 10^° 10^= 10^' 

i [TV- 
Figure 1. The effect of the two evolutionaiy models described in sec- 
tion l2.1l produced by integration of the equations ^6) and 0. (a) Constant 
^1 and invariant Coulomb logarithm, (b) Constant §i and A'^-dependent 
Coulomb logarithm. The multiple lines denote clusters with A^o determined 
by successive factors of four, between 1024 (labeled Ik) and 65536 (labeled 
64k). The different points at which mass loss begins are a physical effect of 
the variation of tf[, with A'^, as for a Plummer Model the core collapse time 
is approximated by a constant number (for this illustration, 20) of initial (.h . 
For the illustrations, f = 0.1 and gi = 0.01. 



these stars will have positive kinetic energy, obtained by accelera- 
tion during the close encounters of two or more bodies in the core 
( lHenonlll96(]h . The probability of an encounter leading to the es- 
cape of a star will be a function of number density and velocity 
dispersion of the stars, with most escapers predicted to come from 
high density regions (namely, the central core). Thus, it is logi- 
cal to suggest that an increased core density will inevitably lead 
to incr eased escape rate. However, the results of ISaumgardt et al.l 
( 120020 show this scaling of ^ is small compared to its overall value, 
so as a first approximation we choose to model the isolated mass 
loss via equation ^ with constant, best-fitting ^i for clusters of 
any iVo . 



2.2 Clusters In a Steady Tidal Field 

If a gravitationally bound object is immersed in an external tidal 
field (i.e. that experienced by a globular cluster orbiting a galaxy). 



every star within the cluster will experience an acceleration from 
the galactic potential. However, stars whose orbits take them signif- 
icantly closer to or further from the central galaxy will experience 
an increasingly different acceleration rel ative to the mean ac celera- 
tion over the cluster. Thus, suggested bv lvon Hoemeii(ll957l) . there 
exists a distance whereupon the tidal acceleration due to the galac- 
tic potential acting on a star exceeds the acceleration between star 
and cluster centre. At this distance, a star would become unbound 
from t he clu ster, and begin to orbit the parent galaxy independently. 
I King 119620 showed that this limiting distance defines the radius of 
the Jacobi surface, which is for a point-mass galactic potential and 
including the centrifugal force given by 



-Rg 



V 3A/g ) 



(18) 



The dynamical effect of the Jacobi radius will be the addition 
of an upper limit to the size to which a cluster can expand. By ap- 
plying our argument that the flux of energy is constant through any 
surface, the perpetuation of this flux requires stars carrying energy 
to cross the Jacobi radius, and hence exit the cluster The result will 
be accelerated mass-loss, significantly reducing the lifetime before 
final dissolution as compared to an isolated cluster. We note also 
that, on account of the A^^' "^ dependence in equation dlSI l, escap- 
ing stars will cause a reduction in the Jacobi radius, which (for 
self-similar evolution) will cause a shrinking of the entire cluster 
(i.e. reduction of all Lagrange radii). 

In a tidal field we require a new definition of ^ to account for 
th e loss of stars t ransiting the Jacobi radius. Following the method 
of lHenon ( 1196 ID , we assume that the existence of such a bound- 
ary will result in a reduced escape velocity, as stars must no longer 
reach infinity. Using the definition of the escape velocity of a tidally 
limited cluster (ft idai) and the e scape velocity from an isolated 
cluster (Diso) from ISpitzeiJ ( Il987h . we find the ratio of escape ve- 
locities 72. ex iityai/^L P C r/rj. 

From this ratio, iGieles & Baumgardd ( 120081) integrated 
Maxwellian velocity distributions to find the fraction of stars with 
sufficient velocity to escape for different ratios 7?.. They showed a 
good fit was given by ^ oc exp (107?.), although for further analysis 
we take a simplified form. 



C = 7C 



2L 

Til 



(19) 



where TZi defines a reference value of TZ and 2 is a constant power. 
We obtain the factor of (3/5) (" present in equation ( I19t as a result of 
Henon's interpretation of evolution occurring with constant mean 
density within clusters. The logarithmic slope relating 7?. to TV is 
given by, 



dln(7?.) dlnr dlnrj 



din TV 



din TV 



(20) 
(21) 



din TV 

3 i 

where we have eliminated d In r/d In N with equation (|9) and de- 
rived d\nri/d\nN from equation dlSt . From equation (I21t it is 
evident that for ^ = (3/5)(^, dln(7?.)/dln A'^ = and so r oc n. 
It is thus apparent (assuming TV < 0) that if ^ > (3/5)^, r will 
shrink faster than rj. Likewise, if ^ < (3/5) (", r will shrink slower 
than rj, with the result that (3/5) ^ is a critical rate whereupon the 
cluster shrinks with constant density. 

We leave the value of z in equation \\9\ to be determined by 
fitting, although note it's value will effect the variation of ^ with 
TZ. The most visible evolutionary effect of the value of z in this 
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regime is to vary the beliaviour of 7?. against A'^, tlie consequences 
and effects of whicli are explained in detail in Appendix IaI 



2.3 Interpretation of the Jacob! Surface 

Tlius far, we have considered the escape of individual stars to be a 
function solely of energy (that is to say, stars with sufficient velocity 
will always escape). It is however possible that stars with energy 
slightly in excess of the critical energy remain bound, owing to 
geometric constraints. Specifically, escape is only possible through 
'apertures' around the Li and L2 Lagrange points, over which the 
gradient of the potential is sufficiently shallow for stars with energy 
only slightly greater than critical to escape. It follows that stars with 
greater energy can escape over a larger region, whilst those with 

less energy will experience a smaller escape aperture. 

As a result of these apertures, iFukushige & Heggid ( I2OO0I) 
characterised the relationship between tesc (the time taken for 
any given star to escape) and the excess energy {E — -Edit) 
(where _Ecrit is the exact energy r equired for e s cape) to be icsc oc 
[-Ecrit/(-E- J5cnt)]^, from which iBaumgardtl ( I2OOII) obtained a 
form for the lifetime of a cluster. 



the variation of ^i (the mass loss of isolated clusters) will be (ap- 
proximately) opposite to that of ^tidai, and thus let (1 — P{TZ, N)) 
represent the 7?. and A'' dependence in ^1. We therefore write. 



(22) 
(23) 



Lev CX t|-h^esc 5 

f N[n{'yNi) Y~'^ 

where A^i is a separate scaling value of N that shall be determined 
from comparisons to A'^-body simulations. It follows that, since a 
star with E > Ea-it will experience a delay before finally exceed- 
ing the Jacobi radius, we should account for these stragglers in our 
evaporation rate. We therefore combine the escape timescale for 
stars (equationl23t. with the rate through which evaporation causes 
stars to evaporate from the cluster (equation [19), and hence obtain 



C = -CP(7^,iV), 

where 



(24) 



(25) 



an improved form for the dimensionless evaporation rate of a clus- 
ter in a tidal field. 



3 A FULL MODEL FOR THE LIFE CYCLE 

Sections [2. 1 | and |23] established physical arguments to predict evo- 
lution of the mass and radius of clusters undergoing either expan- 
sion or contraction. It is apparent that a realistic cluster starting it's 
evolution with r <^ ri will expand to fill it's lacobi radius, eventu- 
ally becoming tidally limited and thus experiencing both regimes. 
We therefore attempt to merge these two formulae into a single uni- 
fied model, encompassing the entire lifetime of star clusters. 

For clarity, we shall refer to the tidally limited escape rate 
(equation I24l> as ^tidai. We also assume that both mechanisms for 
mass loss are viable channels through which stars can escape over 
the entire lifespan of the cluster. However, the extent to which the 
two mechanisms are significant will vary as the cluster expands, 
since early times are likely to be dominated by tidal field indepen- 
dent mass-loss mechanisms and late times by tidal mechanisms. To 
account for this variation we take the factor P{TZ, N) from ^tidai 
(used to represent the 7?. and A'^ dependencies in this term), and 
which varies such that < P{TZ, N) < 1. We now assume that 



^ = ei(i-^(KA^)) + c..dai 

= Cl(l-p(7^,7V)) + |cn7^,Af). 



(26) 

(27) 



Although equation ( I27t represents an essentially simple interpreta- 
tion of the major mass-loss mechanisms present in star clusters, we 
believe an interpretation of this nature is a sufficient description if 
the two mechanisms are both present throughout the lifetime. 

Using equation J27t . we now outline the operation of our nu- 
merical integration code EMACSS. The principle properties, N 
and r, are recovered by application of a 4th order Runge-Kutta nu- 
merical integration kernel to our defining equations (|6) and (O with 
^ and /i calculated appropriately at each integration step, and con- 
stant (". The duration of each time step is set to be O.ltrh, which we 
find to be a reasonable compromise between speed and accuracy 
in the model. This fraction is appropriate as ^ and ^ are defined 
as instantaneous values for the rate of change of A'^ and r per half- 
mass relaxation time. Thus trh forms an upper limit to time step, 
while an overly short time step is liable to be affected by numerical 
inaccuracies. 

At each integration time ti, our procedure is as follows: 

(i) Characteristic properties of the cluster (rj and irh, required 
for the dynamical evolution), are evaluated from N{ti) and r{ti). 

(ii) ^ (equation |27b is evaluated for instantaneous values of 
N{ti) and r{ti), and used to calculate fj, via equation ^. 

(iii) A Runge-Kutta integration step is applied to equations ^ 
and l|7j using parameters derived in stages (i) and (ii), with these 
properties re-evaluated as appropriate. In this stage N{ti+i) and 
r(ii-fi) are recovered. 

(iv) N{ti+i),r{ti+i), and other properties (as required) are out- 
put. 

(v) Steps (i) through (iv) are repeated until A'^ ^ 200, at which 
point the half-mass crossing time tcr ~ irh, and balanced evolution 
is no longer a valid assumption. 

For simplicity, we introduce TV-body units ( such that G = Nrh — 
—AE = r- = 1: iHeggie & Mathiei3ll98d) , although conversion 
to physical units is trivial (EMACSS, by default, outputs both). 
We use the above procedure to illustrate sections 12.11 and 12.31 In 
figure [2] we demonstrate the evolution of a cluster with ^ defined 
by equations l|19lland l l24t . Meanwhile, figure [3] demonstrates the 
evolution of the system using a ^ defined by equation l |27t , for a 
variety of TVq. 

The form of equation Ml\ , and hence the model EMACSS 
are dependent upon several free parameters (see table [T]- However, 
degeneracy amongst these parameters is such that we can adopt 
an appropriate value for TZi, and interpret the scaling factor A'^i 
to define an 'ideal' clusters for which this value is exactly cor- 
rect. We therefore choose TZi — 0.145 (Henon 1965). We addi- 
tionally note that x, z and A'^i all demonstrate a high degree of 
covariance (see Appendix [At. although are not totally degenerate. 
However, to the level of accuracy for which this model is intended, 
this covariance is sufficient that we can eliminate a further param- 
eter without loss of generality. To this end, we choose x — 0.75 
B aumgardt ( 200 1|), and cho ose a value of 7 = 0.11 from the results 
of iGiersz & Heggid ( 119941) for the Coulomb argument, leaving our 
model dependent only upon ^i,Ni,z and ^. These remaining (free) 
parameters are determined by calibration of our model against a 
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Figure 2. Predicted evolution of star cluster parameters according to the models outlined in section lTzl and computed using the technique outlined in section|4] 
(a) shows the evolution of parameters ^ and ^l in time, scaled to C,, while (b), (c), and (d) show the number of remaining stars, half-mass radius, and ratio of Ti. 
respectively. Early divergence between the models is principally an effect of the scaling with N^~^, whist the late divergence between the models (especially 
noticeable in panel (d)) is an effect of the Coulomb logarithm, which becomes dominant at low N. For the illustrations, f = 0.1, A'^i = 10000, TZi = 0. 145, 
7 = 0.11 and x = 0.75, for A'^ = 65536 clusters with TZq = 1/100. For such clusters, the lifetime i^v ~ 700000 A'^-body times, with core collapse 
occurring after ^ 12000 A^-body times. If these clusters were to consist of O.5Af0 stars, an initial r of Ipc would imply an expected lifetime of ~ 60Gyr, 
with core collapse occuning at IGyr. 



database of A'^-body simulations with differing initial conditions, 
which we describe below. 



4 DESCRIPTION OF THE N-BODY SIMULATIONS 

We produced a number of A'^-body simulations using the coUisional 
fourth order Hermite iV-b ody code Nbody6 JMakino & AarsethI 
ll992l : lAarsethll999ll2003l) . on Intel Core 17 computers. Parallelisa- 
tion and GPU acceleration was provided by NVDIA GTX GeForce 
graphics processing units. Our simulations evolved clusters con- 
taining between A'o = 1024 and A'o = 65536 equal-mass stars, 
with a separate series of simulations carried out at every Interven- 
ing fac tor of two. The clusters were Initially described bv' Plummeil 
( Il91l[) models, with initial rt, = 0.78, and henceforth allowed to 
evolve until A'^ ^ 200 (approximately the time of final dissolution, 

fev). 

We defined the ambient tidal field for each simulation in terms 
of initial ratio of half-mass to Jacobl radius 7?.o, using initial ra- 
tios 1/30 and 1/100. The tidal field Implemented was that of a 
point-mass galaxy, with appropriate Ma and Rg to give the re- 
quired rj via equation ( 1181 ). The numbers of our simulations are 
summarised into table|2] along with several simulat ions of isolated 
clusters kindly provided bv lBaumgardt et al.l ( l2002h . 



Table 1. Definitions of the parameters used to describe cluster evolution in 
EMACSS. 



Variable 



Afi 



Definition 



(^ Fractional change in energy per (half-mass) 

relaxation time. 
5 Dimensionless escape rate. 

fi Dimensionless expansion rate. 

TZ Ratio of radius to Jacobi radius. 

4o Time taken for balanced evolution to begin 

in isolated clusters. 



Parameter Definition 



Dimensionless escape rate of an isolated 

cluster. 

Ratio of radius to Jacobi radius for iHenonI 
h965l) models (7^l = 0.145). 
'Ideal cluster' for which 7?.i is exactly cor- 
rect. 

Scaling of g with N around 7?.i . 
Variation of g with tj-j^ . 
Argument of the coulomb logarithm. 
Time in which core collapse is fully com- 
pleted (> io)- 
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Figure 3. The predicted evolution of star-clusters containing varying A^o, using equation j27t . (a) shows the evolution of parameters g (solid) and /i (dashed) 
in time, scaled to f , while (b), (c), and (d) show the number of remaining stars, half-mass radius, and ratio Ti, respectively. The multiple lines denote clusters 
with A^o determined by successive factors of four, between 1024 (labeled Ik) and 65536 (labeled 64k). The two post-collapse regimes are most visible in (b), 
with early mass loss on account of ^i and late due to ^tidal- A maximum size is achieved at 4/tev ~ 0.45 (dependent on initial 7?.), whereupon /i becomes 
negative and the cluster begins to shrink on account of the decrease in Jacobi radius due to decreasing cluster mass. The final changes in the cluster (after 
around i/icv = 0.9) are likely accounted for by the model breaking down, as here t„ K. trh, and thus C, non-constant. The similarity of evolution in panel (b) 
is likely an effect of weakness of the scaling of g with A'^. These illustrations use the same parameters as figures [T]and|2] 



Table 2. The number of simulations carried out f or each set of condition s 
A'^, 7^0- The isolated simulations were provided bv lBaumgardt et al.|j2002l) . 







7^o 




N 


1/30 


1/100 


Isolated 


1024 


64 


64 


3 


2048 


32 


32 


1 


4096 


16 


16 


1 


8192 


8 


8 


1 


16384 


4 


4 


- 


32768 


2 


2 


- 


65336 


1 


1 


- 



For isolated clusters, unbound stars were defined as those with 
energy in excess of that required for escape (i.e. positive energy). 
These escaping stars remained within the simulation until they 
reached 20r, in order to retain their gravitational influence upon 
bound stars. Meanwhile, in tidally limited clusters, the Jacobi ra- 
dius was calculated iteratively using equation dlSt , and unbound 
stars were defined as those outside rj. These stars were removed 
from the simulation upon exceeding 2r-]. This criterion for escape 
will consider stars with E > E^w but r < rj to be bounqj. In both 

* An alternative, energy based formalism has been presented by 



situations, a list of (bound) particles ordered by radial distance was 
used to determine the half-mass radius. 

We finally measured the energy budget of bound stars at each 
time-step. For this purpose, we computed the 'external' energy (ki- 
netic and potential components of single stars and the centres of 
mass for multiples, iGiersz & Heggiell997h separately from the 'in- 
ternal' energy of particles (that stored in binaries and multiples). 
In addition, we recorded the cumulative energy loss from the clus- 
ter for single star escapers and escaping multiple systems. Using 
these data, we were able to track the change of energy, and hence 
(given our premise of balanced evolution) the energy flux through 
any given radius. 



5 RESULTS 

The procedure for the calibration and verification of our prescrip- 
tion against A'^-body data is divided into three principle sections. 
As a first stage, we begin by examining the flow of energy 



iLee & OstrikeJ l ll987l) . although since we seek to include the effect of es- 
caping stars into our model, we choose instead to apply a radial criteria for 
escape. 
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Figure 4. Energy budget for a representative A^ = 65536 run evolving in a 
tidal field with initial r/rj = 1/100. The black (solid) line corresponds to 
the external energy of stars in the cluster, the red (dotted) line is the internal 
energy, the green (dot-dash) line the cumulative energy of singular escapers 
and the blue (dashed) line the cumulative energy of escaped multiples. The 
early (invariant) relationship between energy and time is the pre-core col- 
lapse evolution, and thus before the start of the energy evolution of the clus- 
ter. The downward 'spikes' in the intemal energy correspond to three body 
interactions, in which a single star is ejected with positive (kinetic) energy. 
Meanwhile, the binary involved in such interactions hardens, before finally 
becoming sufficiently hard to escape, carrying away it's (negative) binding 
energy. Effects of this type can be seen to occur intermittently throughout 
the life cycle. 



throughout the cluster. To this end, in section 15.11 we use the en- 
ergy budget of bound stars to determine E/ E by numerical dif- 
ferentiation of E throughout the evolution. We then calculate i,h 
through equation (|4j, and hence determine (, through rearranging 
equation ([3}. This quantity is the 'driver' defining the evolution of 
a cluster, and hence forms a significant limit on the speed through 
which a cluster evolves and the extent of the total lifetime. 

Section [J!2l contains the results of our calibration of EM ACS S 
against A''-body simulations of isolated clusters. In these cases 
^ = ^1 since rj is infinite and ^tidai = 0. Accordingly, we re- 
cover a value for ^i that best represents the evolution of isolated 
clusters. We then apply EMACSS to our TV-body simulations of 
tidally limited clusters, in section [53] and thus assign values for 
the remaining parameters, A^i and z. Using these values, our pre- 
scription is able to concisely express the evolution of A'^ and r for 
clusters throughout their entire life cycle. 

In each case, our best-fitting parameters are recovered by 
a cus tom designed Markov Chain Monte Carlo (^Metrop olis et al.l 
[l953J) routine. For simplicity, we take a series of sainple data en- 
compassing the entire evolution of our A'^-body simulations, and 
apply a Gaussian likelihood function to compare the samples' 
N and r values to the equivalent TV and r values predicted by 
EMACSS. Uncertainties throughout are calculated from the stan- 
dard deviation of our TV-body simulations, but limited to a mini- 
mum of 10% of the absolute value (the level of accuracy to which 
we believe our model to be valid). By marginalising the parameters 
we obtain best fitting parameters for each quantity in the model, 
and can estimate the uncertainty of each. 



5.1 The flow of energy 

Figure |4] shows the evolution of the three principle components of 
the energy budget typical for our TV-body simulations. The quan- 
tities expressed are the total external energy i5cxt of bound stars 
within the system (kinetic and potential components), the total in- 
temal energy E^m (that in binaries and multiples) and the cumula- 
tive total energy carried away by singular escaping stars (iSscsc) or 
multiples (i?mesc). As the cluster is formed in a nearly isolated state, 
early escapers have positive kinetic energy, since they are ejected 
via close encounters with binary or higher order multiples. Accord- 
ingly, corresponding to each escaping (single) star is a decrease in 
the internal energy of the cluster, as the multiple system involved in 
such encounters will become more tightly bound. Through interac- 
tions, a multiple will eventually be in turn ejected from the cluster, 
and hence carry away it's negative (binding) energy. 

Stars ejected late in the life cycle are typically less energetic 
than those ejected earlier This would arise on two accounts; firstly, 
as demonstrated in figure |4] much energy is lost during expansion, 
and hence the stars of a tidally limited cluster retain less energy. 
Secondly, the escape velocity of a Roche lobe filling cluster is re- 
duced by presence of a Jacobi radius, with the result that outlying 
stars with negative energy can exit the cluster 

We begin our analysis by finding the value of C^ (see section|2]l, 
as this is the defining characteristics upon which our model is built. 
For this purpose, we use the variation of the external energy iSexi as 
a function of time, and sequentially bin these data logarithmically 
in energy. In each bin, we approximate the change in the log energy 
over the log time extent of the bin to be linear, and hence determine 
the gradient d log(— i?)/d log t within each bin. 

We then estimate the mean trh corresponding to each bin via 
equation ©.We have previously chosen to work in terms of TV and 
virial radius r, although note trh is defined in terms of rh. We there- 
fore first make the assumption that rh = r to determine a value for 
trh, before rearranging equation (|3) such that 










dlog(-£0 



d log ti 



(28) 



for an arbitrary bin i. We finally substitute trh, ii, and 
dlog(— i5i)/dlogii into equation l |28t to obtain a distinct value 
of C, corresponding to each bin. For completeness, we also measure 
the value of (" for a relaxation time recovered using the measured 
r\ of TV-body simulations, and compare the resultant two values of 

C- 

Figure|5]shows the variation in C, with time, for each of our ini- 
tial tidal field strengths. Although in each case numerical differenti- 
ation has lead to significant scatter, we find that both measurements 
of C, vary around constant values, (" « 0.111 if trh is measured 
from r and C, ~ 0.105 if ti-h is measured directly from rh. These 
results are si milar to previously calcu lated values for Henon's mod - 
els (C = 0.1 ( lGoodman&Hu3ll989l) . C = 0.14 JGao et al.|[l99ll) . 
C, = 0.0 926 for isolated c lusters, C = 0.0743 for tidally limited 
clusters JGieles et al.ll201 11) ). By contrast, there is some indication 
that the smooth flow does not occur during two stages of evolution; 
the first of these corresponds our earliest (furthest left) data, which 
are significantly lower than the average flow. This is the result of 
pre-collapse evolution, during which energy does not change and 
hence C = 0. The second (toward the end of evolution) corresponds 
to the final break down of the smooth energy flow, as at these times 
trh ~ ici and the redistribution of energy becomes increasingly ran- 
dom throughout star clusters. 
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Figure 5. Variation in (^ as a function of time for a vaiiety of rans with different A^o- Only mns with A^o ^ 8192 stars are plotted, as the stochastic nature of 
the energy makes the scatter of smaller runs progressively worse. Filled shapes and the dashed line correspond to measurements of (^ where we have assumed 
rjj = r, while unfilled shapes (and the solid line) corresponds to measurements of ( where r^ is measured directly from A^-body simulations. Panel (a) shows 
the runs carried out with an initial TZ = 1/100 while panel (b) has initial TZ = 1/30. The early (systematically low) points are likely caused by no energy 
loss occurring during core collapse. There is a possible downward slope at late times, which may be a result of the break down of balanced evolution. 



5.2 Evolution of isolated clusters 

We now begin our investigation of EMACSS's ability to reproduce 
the evolution of star clusters, beginning with the simple case of a 
cluster evolving in isolation. In the previous section we demon- 
strated (" — O.lOqj, which we now adopt as the energy flow driv- 
ing the evolution of A'^ and r. We use the numerical integrator de- 
scribed in section|3]to evolve our model cluster, defining the galaxy 
around which our simulated cluster is 'orbiting' to have no mass, 
and hence our cluster to be isolated from tidal effects. The resul- 
tant evolutionary tracks of A'^ and r are compared against those of 
A''-body simulations of isolated clusters. 

For calibrating our prescription we use ^i as a free parame- 
ter, and introduce two further parameters, /m and fr, the fractional 
change in A^ and r during core-collapse. Th e time taken for this 
collapse (assuming an initial IPlummeg ( Il91in sp here) has been re- 
ported to be 17.6 i nitial relaxation times (trh.o) JTakahashil 1 1 9951 : 
IPrukier et al.lll999l) , in which we expect relatively small changes 
in A'^ and r. However, in order that the entire collapse (and core 
bounce) is complete when we begin to apply our prescription, we 
choose to use a time offset of 20iih,o, and define /m and fr to be 
the fractional changes that have occurred within the cluster prior to 
this time. Naturally, it is therefore also possible to project our pre- 
scription backward until the time of core collapse itself, and hence 
model the entire post-collapse evolution of a star cluster. We find 
values for /m and fr from linear interpolation of Baumgardt's A'^- 
body data at 20f,h,o, and calibrate our prescription using the Monte 
Carlo algorithm described in section |3] 

Figure |6] shows the evolutionary tracks predicted by our best- 
fitting parameters for the entire lifetime of the cluster Our best fit- 
ting models occur where ^i = 0.0141, consistent previous fin dings 
(Ci « O.Ol.'Heggie & A arsetlj[l992l : iBau mgardt et al.l l2002l) . Our 
interpolation demonstrates a mean (A'^ independent) expansion dur- 
ing core-collapse fr — 1.81, with the majority of this expansion 
occurring immediately after the predicted 17.6frh,o core collapse 



Using (^ measured from A'^-body simulations. 



time. Meanwhile, mass-loss occurs such that /at = 0.95. From fig- 
ure |6l it is apparent that this mass-loss begins gradually, building 
up over a period of time to approximately constant ^i observed for 
the majority of the life cycle. It is clear therefore that backward 
projection of our prescription will result in an instantaneous start to 
mass-loss, an effect that we do not believe to be physically likely. 
We therefore conclude that the early evolution - directly after core 
collapse - will be poorly described by this prescription 

The principle discrepancies between our best fitting prescrip- 
tion and comparison A'^-body data consist of some mild variation 
of log-gradient d log r/d log A''; while we represent the evaporation 
rate as being constant, the A^-body data disp l ays a more complex 
evolution, a point noted bv lBaumgardt et alj ( 120021) . In this paper 
the authors demonstrate that expansion is not strictly self-similar 
for an isolated cluster, although they find similar scaling between 
expansion and mass-loss at an inner Lagrange radius (10% — 20%) . 
We consider this to be a consequence of an increased rate of ejec- 
tion causing encounters on account of increased core density for 
larger A' systems, and hence interpret the discrepancies in our gra- 
dient to be a result of clusters not being self-similar throughout 
their evolution, although note our current self-similar description 
does provide a good first order interpretation. 

We demonstrate the veracity of our best fitting parameters 
by considering our first approximation of constant Coulomb log- 
arithm. Here, A^ should go as a negative power-law with time 
(equation II 61) ■ with the power i/ defined by equation ( I12t . Using 
( — 0.105 and ^i = 0.0141, we compare the init ial mass-loss rate 
from our results to the mass-loss rates found by iBaumgardt et al.l 
(2002) for our range of A'o. Our measurements give that v = 0.130, 
independent of TVo. This value of u demonstrates excellent agree- 
ment with Baumgardt's values of i/ calculated by the evolution of 
AT (0.13 s; f ^ 0.14 for 1024 s^ A^ ^ 8192), although less 
good agreement with Baumgardt's f calculated by the evolution 
of r (0.02 ^ !/ ^ 0.09 for 1024 ^ A^ sj 8192). We conclude 
therefore that equation ( I16t can be used to calculate the relation- 
ship between A' and time, but that this scaling cannot simultane- 
ously reproduce the evolution of r with comparable precision, as 
TV and r do not scale with the same power-law. 
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Figure 6. (a) The relationship between A'^ and r for clusters of different A^o- (b) The evolution of A'^ (scaled to the initial cluster size A^o)- The scattered 
points denote the evolution of clusters with A^o determined by successive factors of two, between 1024 (labeled Ik) and 8192 (labeled 8k), while the black 
lines are produced using our prescription with our best fitting parameters. Our (best-fitting) pres cription provides a gene ral indication of the rate at which an 
isolated cluster expands and evaporates, although does not reproduce the crossing point found bv lBaumgardt et al.l 112002 1. There is some evidence (explained 
below) that the systematic discrepancies between our prescription and the Af-body data are effects of isolated clusters not being fully self-similar throughout 
the entirety of their evolution. 



We finally examine our definition of to (eguation llSb . which 
we have described as being the time required for a model of non- 
zero ro to join the track predicted by balanced evolution (where 
t = r = Q). Using the value of ^ and C, above, we find to ~ 9.2frh,o, 
a time shorter than the core-collapse of a Plummer profile. This 
value does demonstrate the dependence of the length of the core- 
collapse upon initial conditions, but does not give any insight into 
the subsequent global evolution. 

5.3 Evolution in a tidal field 

A cluster in a tidal field will first experience an expansion phase, 
before becoming significantly affected by the tidal field and en- 
tering a contraction phase. We assume the expansion phase to be 
similar to that which we have previously demonstrated for an iso- 
lated cluster, and hence expect a value of ^i ~ 0.014. We initially 
set (^ = 0.111 (using our assumption that rh = r when calculating 
the relaxation time) and fee = 20.0fih,o- 

For tidally limited clusters, we expect the mass-loss prior to 
core collapse to be more significant than that found for isolated 
clusters. From the evidence of equation J24t , we expect this mass- 
loss to be a function of TZ, and therefore allow a unique /jv and 
fr for each simulation. Accordingly, we measure these values of 
/jv and fr by linear interpolation of our A'^-body data at 20irh,o> 
and consequently recover the specific fractional changes in each 
series of simulated clusters during core-collapse. These fractional 
changes in the size and mass of clusters are demonstrated in fig- 
ure |7] 

The extent of mass-loss and expansion during core-collapse 
shown in figure |7] reveals a marginally greater expansion but 
lower mass-loss for higher TZo- In each case, simulated clusters in 
stronger initial tidal fields expand less, but loose more mass than 
those in weaker tidal fields. Additionally, figure |7] suggests Jm is 
only weakly dependant on A^o, implying that the fraction of stars 
ejected prior to core collapse is not significantly dependant upon 
the total number. We use these direct measurements of expansion 
and mass-loss for the purposes of our fitting. 
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Figure 7. Mass-loss (/jv) and expansion (/r) of clusters pre core-collapse. 
The most significant dependency is that with TZo (on account of the overall 
escape energy). There is a small A^q dependency visible in /jy, although no 
particular relationship between fr and A^q. 



We employ our Monte Carlo algorithm to simultaneously fit 
our prescription to each of our tidally limited A-body simulations 
(table |2j, comparing the evolution of ^i, A and r. The resultant 
joint posterior probability density demonstrates a narrow peak (as 
shown in Appendix|Cj. This is most probably a result of covariance 
within our data, since each datum will be highly correlated to those 
around it. We therefore take our repeated database of TV = 4096 
and N = 8192 simulations, and independently fit each of these 
with our prescription, producing a unique value of A^i and z. We 
find the standard deviation of these individual A^i fits to be about 
40% of the mean value, while the standard deviation of the z fits 
is about 15%. We hence interpret these standard deviations as esti- 
mates of the uncertainty on these parameters. 

The best fitting parameters we recover are shown in table |3] 
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whilst we over-plot the predicted evolution from our prescription 
to A'^-body data in figure[8](dotted lines). For the sake of complete- 
ness, we additionally show example joint and marginal posteriors 
distributions recovered by our Markov code in Appendix |C] and 
additional compari sons against the simpler analytic models from 
iGieles et alj J201lh in Appendix H) 

The best-fitting values we recover are similar to those previous 
predicted. We find that for the tidally limited cluster, ^i = 0.0142, 
similar to that required by the isolated cluster, and A'^i = 38252, 
comparable to, although smaller than, that previously anticipated 
(A'^i ~ 10^, although d erived from models containing a mass func- 
tion lGieles et alj201 ih . We additionally find z = 1.61 (remarkably 
close to z = 1.5, the value chosen for convenience in the same 
paper). There is evidence of some covariance between A'^i and z, 
which probably arises due to the manner through which these two 
parameters define the total lifetime of clusters. The total lifetime 
(see Appendix [At is a function of z, Ni, and x such that the rela- 
tionship between z and TVi plays a substantial part in the average 
rate of mass-loss throughout the lifetime. Accordingly, a variation 
in one is likely to cause a comparable variation in another, which 
will be expressed as some degeneracy. 

We find that our model is able to reproduce the evolution of A*' 
and r, and the correct behaviour in the N — TZ plane. Figures[8la) 
and (b) show that the early expansion phase, transition phase, and 
evolution in the tidally limited regime are all well reproduced by 
our prescription. However, some (small) systematic offsets are vis- 
ible in our prescription - typically observed around the point of 
transition at which the tidal field becomes significant. This is the 
result of two effects: first the none self-similar evolution discussed 
in section [5^ and secondly the simplistic manner in which ^tidai 
and ^1 are combined in equation ( |27t . Nonetheless, this offset does 
not degrade the accuracy of this model by a significant (> 10%) 
margin, and is hence acceptable for our purpose. 

It is interesting to note that in both figures|8la) and (b) TZ{N) 
is not constant with decreasing A'^ in the asymptotic regime (the 
line upon which all our evolutionary tracks merge is not horizon- 
tal). This suggests that in this stage of evolution the density of the 
cluster is decreasing as stars escape. The cause of this effect is 
demonstrated in appendix lAl - an A'^ dependent ^ for a; < 1. We 
find however that such an effect is not explicitly demonstrated in 
the evolution of rn/rj (instead, we find once again approximately 
constant density), although we leave the investigation of this phe- 
nomenon to a subsequent paper. 

Panels (c) and (d) show the time-evolution of A'^ for our clus- 
ters. In both cases, the decrease in TV is adequately reproduced, 
although some particular series of simulations show systematic off- 
sets. These tend to arise on account of small variations in the time 
taken for core collapse, and are therefore likely to be a statistical 
effect of our A^-body simulations. Similar effects are also observed 
in the variation of r (panels e and f), once again most likely on 
account of random perturbations in the initial distribution of stars. 

Our prescription is generally successful in producing our A'^- 
body data to within the standard deviation of our series of A'^-body 
simulations (see figure |9]l. Despite this, our TV-body data shows 
some overlaid noise (i.e. oscillatory behaviour), that is particularly 
apparent for the evolution of r when compared to that of TV. This 
noise becomes especially significant in later evolution (e.g. contrac- 
tion), whilst expansion is comparatively smooth. We believe this 
effect is indicative of the ejection of stars, since every ejection of 
a single star or binary will be accompanied by a change (jump) in 
the total energy of the cluster. Accordingly, because the total en- 
ergy is lower in later stages of the evolution, these jumps are more 
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Figure 9. The median values (solid lines) and standard deviations (shaded 
area) of our repeated simulations of A^o =1024, 4096, and 16384 clusters. 
The standard deviation is substantially greater for smaller A^o simulation, as 
statistical noise is more significant when the number of stars is small. Ac- 
cordingly, the standard deviation of our simulations increases with t (and 
hence decreases with increasing A'^). We therefore require that the prescrip- 
tion describes higher A^o simulations well, while the significance of our 
results at lower iVo is reduced. In order that undue significance is not given 
to the early evolution (with high TV), our limiting 10% error is especially 
applicable in these regimes. The smooth (black) lines over-plotted are our 
individual best fits to these simulations. 



significant. These events demonstrate no particular TV dependency, 
and, as an effect corresponding to randomly occurring events, are 
not naturally reproduced by EMACSS. 



6 SUMMARY AND CONCLUSIONS 

We have developed a prescription and code to reproduce the post 
core-collapse evol ution of t he ma s s and radius of star clusters, built 
upon the w ork of HenonI 1 196lL 119650 , whose combination was 
proposed in lGieles et al.l (1201 ih . The prescription is fundamentally 
simple, considering a cluster of equal-mass stars in the tidal field of 
a point-mass galaxy. 
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Figure 8. The evolution of star clusters in point mass tidal fields. The scattered points denote the evolution of clusters with A^o determined by successive 
factors of two, between 1024 (labeled Ik) and 65536 (labeled 64k), while the black (dotted) lines are produced using our prescription with our best fitting 
parameters. The left hand column corresponds to an initial TZq of 1/100, while for the right hand T^o = 1/30. (a,b) TV against 7^. The movement in the 
N — TZ plane is consistent for all A^-body simulations; the tracks converge once the cluster is tidally limited, with expansion dominant until this transition. 
(c,d) Af as a function of time. (e,f) Radius r as a function of time. The expansion and contraction stages are accurately reproduced, although the fit for the 
model is less good during the transition between these regimes (i.e. peak r). The oscillations observed in the later evolution of r are possibly indications of 
gravothermal oscillations or stochastic noise in the core. The model once again fails for low N. For 0.5Mq stars, an initial r of Ipc would give an expected 
hfetime of 35Gyrforthe N = 1024, TIq = 1/100 simulations (core collapse after 260Myr), and 55Gyrforthe N = 65536, 7^o = 1/100 simulation (core 
collapse after IGyr). Meanwhile, the expected lifetime of the N = 1024, T^o = 1/30 simulations is around 6Gyr, and that of the N = 65536, TZq = 1/30 
simulation is around llGyr 
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Table 3. Parameters defining the EMACSS code. Those of the left are fixed at literature values, while those on the right are recovered by a Markov Chain 
Monte Carlo technique. The differences in f are present on account of the two different definitions of r used in the expression for half-mass relaxation time. 



Fixed 
7^1 



Markov Chain 

5i Ni 



Isolated 0.11 
Tidal 0.11 0.145 



0.75 



20 
20 



0.105 0.0141 

0.111 0.0142 38252 



1.61 



We allow an approximate form for the mass-loss during a 
cluster's expansion phase, followed by a smooth transition into 
a tidally limited pha s e. In this tidal phase we use the method of 
iGieles & Baumgardtl ( 120081) to calculate the tidally-induced mass- 
loss. We additionally include the variation of the Coulomb loga- 
rithm, which has particular influence in late evolution. As a con- 
sequence, our model encompasses both expansion and contraction 
phases, and is applicable throughout the entire post core-collapse 
evolution of star clusters. Thus, we have unified, and expanded 
upon, the models of Henon. 

We rely upon the axiom that the fractional conduction of en- 
ergy from the core is constant per relaxation time. It is demon- 
strated that this is indeed true, for the majority of the cluster's life- 
time, and that star clusters exist for the majority of their lifetime 
in a state of balanced evolution. From our direct measurements, 
we show that (" ~ 0.1, in agreement with the values predicted by 
Henon's models. This measured value is then used as a foundation 
for our prescription, which, we demonstrate, can be used to model 
the evolution of clusters of equal-mass stars both in isolation, and 
in the tidal field of a point-mass galaxy. We have shown that our 
prescription is generally successful, albeit with certain regions in 
which our prescription accuracy experiences systematic offsets. 

For the isolated regime, we find that EMACSS can be made 
to approximately reproduce the simulated evolution of a star clus- 
ter for the majority of the lifetime, although also demonstrates that 
the evolution is not strictly self-similar. We additionally find ef- 
fects associated with a gradual start to mass-loss, and the failure of 
the axiom of constant (^ during the final dissolution of the cluster. 
Meanwhile, for clusters in tidal fields, we find that the EMACSS 
code can successfully describe all stages of evolution, in terms of 
A'^ and r. There is some stochastic variation in the time taken for 
balanced evolution to begin in our TV-body data, although note that 
if this offset is accounted for our prescription is successful. 

In conclusion, we have added a description of variable mass 
loss to the models of Henon, and have developed EVOLVE Me A 
Cluster of StarS, a numerical code that uses our description 
to predict the evolution of a star cluster. We have calibrated this 
code to recover A*' and r to within 10% for the full life cycle of 
equal mass clusters. Although this description does not consider 
any of the complexity arising from full relaxation mechanisms or 
oscillations within the core, these effects appear to be smoothly 
variant over a half-mass relaxation timescale. Thus, our axiom of 
a constant change in energy per irh, and the dynamical mass-loss 
and expansion rates derived from this, are shown to be conducive 
to predicting the evolution of a star cluster. 



6.1 Future Development of our Code 

Our first iteration of EMACSS is publicly availablqj. The code 
predicts the evolution of mass and radius for a virialised cluster 



https://github.com/emacss 



of equal mass stars, throughout it's lifetime. As the dynamical re- 
laxation effects of the cluster are dominant throughout significant 
periods of a realistic cluster's life cycle, we believe this simple code 
will reproduce realistic star clusters to a good degree of accuracy. 

The inaccuracy of our model at very low (A'^ < 200) numbers 
of stars is acceptable, since observation of clusters on this scale is 
unlikely. Moreover, these small clusters are comparatively fast to 
model via other means. This failure does however define a lower 
limit for the accuracy of EMACSS, while our upper limit is less 
certain. 

Despite these successes, our code remains fundamentally sim- 
ple, with many physical effects neglected. Immediately apparent is 
the lack of a mass function; with a realistic mass function, mass 
segregation is likely to increase the speed of the evolution of a 
clusteilj. Furthermore, preferential ejection of particular (low mass) 
stellar types will further complicate the evolution of a cluster by 
changing the mean mass of stars over time. Thus, clusters with a 
high mass range are likely to be worse reproduced than flatter mass 
functions, and will require modifications to our code. 

The parameters produced are limited to mass and a virial ra- 
dius, which we use as an approximate representation of half-mass 
radius. The approximation of half-mass radius is made through the 
assumption that the energy form factor, k, does not vary through- 
out a cluster's life cycle. It follows that by accounting for this 
form factor we would be able to reproduce the half-mass radius, 
although conversion to physical observables - such as half-light ra- 
dius - would require additional consideration of the stellar mass 
function and exact form of density profile. 

We have also considered only simple formulae to describe es- 
cape mechanisms within the cluster (e.g. our single factor form for 
escape rate in the isolated regime). In this way, we have overlooked 
several factors and mechanisms in the escape process, which may 
account for some of the discrepancies noted between our predic- 
tions and A-body simulations. A more detailed study of these 
mechanisms is forthcoming, with the intention of improving the 
fit of our model throughout the complete life cycle. 

We finally note that our prescription interprets the mass-loss 
due to a tidal field as an effect dependent upon a limiting Jacobi 
surface. Once again, progress is forthcoming on a more precise 
realisation of tidal fields in which the globular cluster exists (e.g. 
iLamers et al.ll201 (J). additional potentials and orbital eccentricity. 
It follows that future implementations of ou r code will incorporate 
these more versatile tidal fi eld descriptions {[Tanikawa & Fukushigd 
uOlduRenaud et al.ll201 ih , and hence allow a greater range of situ- 
ations to be studied. In conjunction with realistic (and evolving) 
mass functions, these improvements could significantly improve 
our results for realistic cluster simulations, with a potential exten- 
sion into population studies of the local universe's globular clusters. 



^ The most likely effect is an increase of i^ , since simulations of clusters 
with mass functions are seen to evolve faster uGieles et al.|j2O10l) found that 
clusters with mass functions typical of globular clusters are well described 
whenC Ri 0.2. 
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APPENDIX A: A GENERAL SOLUTION TO 7^ - iV AND 

N{t) FOR ANY X AND Z 

Although the principal objective of this study is to develop a nu- 
merical prescription for the evolution of clusters, we nonetheless 
explore the evolution analytically to recover (general) relationships 
between our fitting parameters. In particular, we explore the re- 
lationship between x and z, as both these powers are present in 
Ctidai- For the sake of a first order approximation we consider the 
Coulomb logarithm and k to be constant, and hence rewrite equa- 
tion l24l as 

We can now apply equation l lAlb to equation l l21b and obtain 



diV 



n 

N 



5 f Tl_ 
3 V^i 



(A2) 



This can be written as Bernoulli's differential equation, taking the 
form 



^ + p(iv)7^ = Q(iv)7^^-^ 

where 

V{N) = -(5/3)/A^ 
Q{N) = aN""^ 



(A3) 

(A4) 
(A5) 



with a = —{5/3)TZ^ NI ^. We now use a variable substitution 
U = W to get 



+ riN)TZ-' = Q{N). 



1 dW 
zdN 
which can be solved with the correct constant of integration 



M = TZoN^ 



(5/3)3 



(A6) 



(A7) 
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Solving equation l lA6t with l lA7t we find 

x-l / r AT -] (,5/3)z-x+l 



n- 






N_ 



7^o 



7^l 






-,(5/3)3\ 1/^ 



wliere A = ([2/3]^ + [2/5][l - j]/[2/3]z) 

etalJjloTTh 



7?.o = we find equation (A8) of lGieles 

7^l 



7^ = 



1 + 4(1 



0) 



2/3 



iVi 



1- 



(A8) 
For z = 3/2 and 

2/3 






-|(7/2)-a, 



(A9) 

Equation l |A8l ) converges in thie tidal regime to the power-law rela- 
tion 

x-l 

(AlO) 



Substituting this into equation iA2t we find that the logarithmic 
slope around TVi is 

dln7^ 5, 



dlniV 3 

X 



(l-A) 
1 



(All) 



(A 12) 



Thus, the assumed values used in the model of lGieles et alj 1 120 llh 
always has a logarithmic slope (2/3)(a; — 1) = —1/6. 
The value of S, converges to 



^~ 5A' 

while the relaxation time varies according to 



(A13) 



Lrh — f-r 



3/2 



(A14) 



Through combination of equations (|6j, ( IA13t and ( IA14t . we hence 
find that 

TV 



N 



< 



irh 



3 iVi 

5 irh.l 



l3/(2z) 



Til 

iVi 



O trh.l 



-3/2 

iV 

ivT 



3(l-^)/(2z) 



(A 15) 
(A16) 
(A17) 



Combinations of (x, z) can be found such that 3(1 — x)/{2z) = 
0.25 (i.e. {x,z) = (1/2,3); (3/4,3/2); (5/6, 1)). These will all 
give the same value for A as well, but smaller values for z lead to 
slightly shorter lifetimes (equation [ATT). This can also be under- 
stood from the asymptotic behaviour in the TZ — N plane: smaller 
z means smaller TZ which means shorter relaxation time. Since all 
combinations of x and z will eventually converge to the same 5, 
this implies the total lifetimes will be shorter. 



APPENDIX B: COMPARISON TO PREVIOUS MODELS 

In order to demonstrate the improvement of our new prescrip- 
tion when compared to previous models, we attempted comparable 
Markov Chain analyses for a series of models of increasing com- 
plexity. Accordingly, each model increases the number of free pa- 
rameters whose values are to be recovered. The models we attempt 
are as follows: 




100000 200000 300000 400000 500000 500000 



t\N- 



Figure Bl. Different iterations of ^ against A^-b ody data. In case (i ), the 
parameters used for the model are those quoted in lGieles et alj J201lh . For 
the remaining cases, the models are generated using the 'optimal' (best- 
fitting) parameters as recovered by a Markov chain. The parameters used 
are given in table IB II while the forms for § are given in appendix IB] 



Table Bl. Parameters used for the 'best-fitti ng models' shown in figure IbTI 
The parameters for case (i) are taken from iGieles et al.l 11201 1 ), while the 
remaining models have parameters recovered by a Markov Chain Monte 
Carlo method. A comparable Markov Chain analysis was attempted for 
model (i), though failed to converge. The Markov Chain used for model 
(ii) also did not converge to a single value, but instead produces several 
'best-fitting' sets of parameters. 



Parameter 


(i) 


(ii) 


(iii) 


(Equation |27) 


c 


0.08 


0.0638 


0.104 


105 


«1 


- 


0.0121 


0.0085 


0.0161 


7 


0.11 


0.11 


0.11 


0.11 


7^l 


0.145 


0.145 


0.145 


0.145 


Ni 


- 


- 


54721 


26689 


X 


- 


- 


0.566 


0.731 


z 


1.5 


1.5 


2.67 


1.75 
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(i) Model of loieles et alj JJoTIh . 

This model c omprises the simpl e description of ^ presented in 
equation (11) of lGieles et alJ ( 120 llh . namely 



^ = 7C 



t-cr 
tcr,l 

5^ It^i 



3/2 



(Bl) 



where icr is the (half-mass) crossing time. This model provides a 
description only of tidally induced mass loss, and is characterised 
by 4 parameters (C, 7?.i, /n and /r). 

(ii) Model of iGieles et alj ( 120110 with constant isolated mass 
loss. 

For this iteration, a second term is included in ^ to account for 
the mass-loss in the isolated regime. Hence, ^ is written as 



^Mi 



+ Ci 



(B2) 



and includes the additional free parameter ^i. We have also re- 
placed the factor of 3/2 with z which we now allow to be free, and 
fix TZi — 0.145 as this parameter is heavily degenerate with the 
combination of z and (". The effects of this variable z (as opposed 
to a variable 7?.i are discussed in appendix lAt. 

(iii) Appendix model of lGieles et al.l ( 120111) with constant iso- 
lated mass loss. 

Thi s form for ^ is that presented in Appendix A of lGieles et al.l 
( 12011 1). but with the inclusion of an additional term to describe the 
isolated mass loss. We hence use a modified form of iGieles et al.l 
( I2OIII) equation (A4), such that, 



^-H^ 



JVlniVi 
NilnN 



+ Ci 



(B3) 



where we have additionally included the In 7A'^ dependence of the 
Coulomb logarithm. We now have free parameters (", ,^0, -'Vi, x, z, 
/n and /r. 

We use our Markov fitting code to fit the above models to a 
single sample of iV-body data (TV = 65536, 7^l = 1/100), and 
over-plot the resultant evolution for the 'optimally' fitting param- 
eters for each iteration of our prescription in figure Ell We also 
plot the form from our full prescription, in which ^ is defined by 
equation l l27t . and give our best-fitting parameters for each case in 
table IB 1 1 We note however that (i), the simplest version did not 
converge, and hence did not produce any useful best-fitting param- 
eters. In addition, model (ii) appears to have no single value upon 
which the model has converged, but instead several values between 
which the Markov chain has oscillated. 

In each case we have allowed all possible parameters to be 
free. Accordingly, degeneracies between certain parameter (e.g. x, 
Ni, z) are clearly visible in the joint posterior probabilities. 



overall fitting to clusters in tidal fields. The remaining figures show 
posteriors recovered for the simpler models described in Appendix 
B; figure C3 shows the posterior for fitting model (ii), figure C4 is 
the posterior from fitting model (iii), and figure C5 is the posterior 
of fitting out full prescription to a single A'^o = 65536, TZ = 1/100 
simulated cluster. 



APPENDIX C: POSTERIOR PROBABILITY 
DISTRIBUTIONS RECOVERED BY MARKOV CHAINS 

The posterior probability densities produced by our Markov chain 
fitting are shown in this appendix, for a variety of situations in 
which we have tested our prescription. As such, the most probable 
values for parameters, and an estimate of uncertainty are present 
in the histograms, whilst degeneracies are visible in the probabil- 
ity density maps. Figure CI is shows the posterior recovered for 
our overall fitting to isolated clusters, while figure C2 shows our 
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0.01410 



0.01411 



Figure CI. Histogram of the posterior probability for the ^i of an isolated cluster The histogram is normalised to a maximum peak value of 1. We believe the 
mode value to be accurate, although the width of the posterior is likely to be underestimated on account of covariance within our data. 



^ 



38254 














m 




38253 






38252 
1.6102 














1.6109 










m 


• 




1.6108 








1.0 














0.8 














0.6 












^ 


0.4 












^ 


0.2 












. 




1 






1 



0.0142255 0.014226538252 



38253 



1.6109 

Z 



Figure C2. Joint and marginal posterior probabilities of ^i, A'l and z for the complete prescription, comparing the evolution of N{t) and r{t). The pre core 
collapse mass-loss and expansion, /jm and /r, are measured from A'^-body data, and tec = 20.0trii.o- The figure is generated from a single Markov chain 
describing all our range of A^o for ^0 = 1/100, 1/30, and exhibits a very narrow peak on account of covariance within our data. 
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Figure C3. Joint and marginal posterior probabilities of C, and ^i for model (ii) in Appendix IbI Parameter fee is assigned the value 20.0trh.o. ^nd f^ and 
fr are measured from A^-body data. The properties compared are the evolution of A'^ and r. The Markov chain has not converged in this plot, and is instead 
wandering between several sets of 'best-fitting' values. 
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Figure C4. Joint and marginal posterior probabilities of C,, ^i,Ni, x, z for model (iii) in Appendix IB] Parameter tec is assigned the value 20.0trti.o> ^nd /jy 
and fr are measured from Af-body data. The properties compared are the evolution of A^ and r. 
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Figure C5. Joint and marginal posterior probabilities of C^, ^i, A^i, x, and z for the complete prescription. Parameter tec is assigned the value 20.0trh^o, and 
/at and /r are measured directly from Af-body data. The properties compared are the evolution of N and r. Degeneracies between x, N\, ^i and z are clearly 
visible. 



